Geospatial Assignment 2
Irache Vergara, Gretchen Moulton, Viktoria Gagua, Daniela Guggenberger
2025-02-04
1 Load libraries
rm(list = ls())
# Library
library(readxl)
library(sf)
library(ggplot2)
library(maps)
library(dplyr)
library(tidyr)
library(readr)
library(rnaturalearth)
library(rnaturalearthdata)
library(osmdata)
library(remotes)
library(afrilearndata)
library(RColorBrewer)
library(gridExtra)2 Load data
# Define the base URL for raw shapefiles on GitHub
github_base <- "/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/data"
# Define file paths for each dataset
pop <- st_read(paste0(github_base, "/ne_110m_populated_places/ne_110m_populated_places.shp"))## Reading layer `ne_110m_populated_places' from data source
## `/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/data/ne_110m_populated_places/ne_110m_populated_places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 243 features and 137 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -175.2206 ymin: -41.29207 xmax: 179.2166 ymax: 64.14346
## Geodetic CRS: WGS 84countries <- st_read(paste0(github_base, "/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp"))## Reading layer `ne_110m_admin_0_countries' from data source
## `/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/data/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 177 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84airports <- st_read(paste0(github_base, "/ne_10m_airports/ne_10m_airports.shp"))## Reading layer `ne_10m_airports' from data source
## `/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/data/ne_10m_airports/ne_10m_airports.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 893 features and 40 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -175.1356 ymin: -53.78147 xmax: 179.1954 ymax: 78.24672
## Geodetic CRS: WGS 84ports <- st_read(paste0(github_base, "/ne_10m_ports/ne_10m_ports.shp"))## Reading layer `ne_10m_ports' from data source
## `/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/data/ne_10m_ports/ne_10m_ports.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1081 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -171.758 ymin: -54.80944 xmax: 179.3094 ymax: 78.22611
## Geodetic CRS: WGS 84locations <- st_read(paste0(github_base, "/ne_10m_populated_places/ne_10m_populated_places.shp"))## Reading layer `ne_10m_populated_places' from data source
## `/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/data/ne_10m_populated_places/ne_10m_populated_places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7342 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -179.59 ymin: -90 xmax: 179.3833 ymax: 82.48332
## Geodetic CRS: WGS 84# Convert and transform spatial data if needed
countries <- st_as_sf(countries)
airports <- st_as_sf(airports)
airports <- st_transform(airports, crs = 4326)
ports <- st_as_sf(ports)
ports <- st_transform(ports, crs = 4326)3 Graphs and Visuals
# Plot the world map
ggplot(data = countries) +
geom_sf(fill = "lightblue", color = "black") + # Fill color and borders
theme_minimal() +
labs(title = "World Map")
#airports on map
ggplot() +
# Plot world map
geom_sf(data = countries, fill = "lightblue", color = "black") +
# Plot airports on top
geom_sf(data = airports, color = "yellow", size = 1) +
theme_minimal() +
labs(title = "World Map with Airports")
#ports on map
ggplot() +
# Plot world map
geom_sf(data = countries, fill = "lightblue", color = "black") +
# Plot airports on top
geom_sf(data = ports, color = "red", size = 1) +
theme_minimal() +
labs(title = "World Map with Ports")
ggplot(data = countries) +
geom_sf(aes(fill = POP_EST), color = "black", size = 0.2) + # Fill based on population
scale_fill_viridis_b(option = "plasma", trans = "log", name = "Population",
breaks = c(1e6, 1e7, 1e8, 1e9), # Define readable breaks
labels = c("1M", "10M", "100M", "1B")) + # Format labels
theme_minimal() +
labs(title = "World Population by Country",
subtitle = "Based on Estimated Population (pop_est)",
caption = "Source: Natural Earth / Your Dataset") +
theme(legend.position = "right", # Move legend to the right side
legend.text = element_text(size = 10), # Increase text size
legend.title = element_text(size = 12, face = "bold"), # Bold legend title
legend.key.height = unit(1, "cm")) 
plot_population_distribution <- function(continent_name) {
continent_data <- countries %>% filter(CONTINENT == continent_name)
ggplot(continent_data, aes(x = reorder(NAME, POP_EST), y = POP_EST)) +
geom_bar(stat = "identity", fill = "steelblue") + # Bar plot for population
labs(title = paste("Population Distribution of", continent_name, "Countries"),
x = "Country",
y = "Population") +
scale_y_continuous(labels = scales::comma) + # Format y-axis with commas
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
}
# List of continents
continents <- c("North America", "South America", "Asia", "Europe", "Africa")
# Apply the function to each continent and create the plots
plot_north_america <- plot_population_distribution("North America")
plot_south_america <- plot_population_distribution("South America")
plot_asia <- plot_population_distribution("Asia")
plot_europe <- plot_population_distribution("Europe")
plot_africa <- plot_population_distribution("Africa")
# Display each plot separately
print(plot_north_america)
print(plot_south_america)
print(plot_asia)
print(plot_europe)
print(plot_africa)
plot_population_distribution <- function(continent_name) {
continent_data <- countries %>%
filter(CONTINENT == continent_name) %>%
mutate(Pop_Percent = (POP_EST / sum(POP_EST, na.rm = TRUE)) * 100) # Calculate % of total
ggplot(continent_data, aes(x = reorder(NAME, Pop_Percent), y = Pop_Percent)) +
geom_bar(stat = "identity", fill = "steelblue") + # Bar plot for % population
labs(title = paste("Population Distribution (%) of", continent_name, "Countries"),
x = "Country",
y = "Percentage of Continent's Population") +
scale_y_continuous(labels = scales::percent_format(scale = 1)) + # Format as percentages
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
}
# List of continents
continents <- c("North America", "South America", "Asia", "Europe", "Africa")
# Apply the function to each continent and create the plots
plot_north_america <- plot_population_distribution("North America")
plot_south_america <- plot_population_distribution("South America")
plot_asia <- plot_population_distribution("Asia")
plot_europe <- plot_population_distribution("Europe")
plot_africa <- plot_population_distribution("Africa")
# Display each plot separately
print(plot_north_america)
print(plot_south_america)
print(plot_asia)
print(plot_europe)
print(plot_africa)
4 Locations
#locations on map
ggplot() +
# Plot world map
geom_sf(data = countries, fill = "lightblue", color = "black") +
# Plot airports on top
geom_sf(data = locations, color = "yellow", size = 1) +
theme_minimal() +
labs(title = "World Map with Locations")
5 Average Distances’ Histogram (Airports)
world <- ne_countries(scale = "medium", returnclass = "sf")
# Spatially join the population table with world to add country info
pop_with_country <- st_join(pop,
world[, c("sovereignt", "continent")],
join = st_intersects)
# Remove any rows that did not join to a country
pop_with_country <- pop_with_country %>%
filter(!is.na(sovereignt))
# For each country, select the top 20 populated locations (using POP_MAX)
top20_pop <- pop_with_country %>%
group_by(sovereignt) %>% # Group by country name
slice_max(order_by = POP_MAX, # Order each group by the POP_MAX column
n = 20, # Keep the top 20
with_ties = FALSE) %>% # Exclude ties beyond the 20th entry
ungroup()5.1 Airports
#assign each airport a country
airports <- st_transform(airports, st_crs(world))
airports_with_country <- st_join(airports, world, join = st_intersects)
airports_with_country$sovereignt## [1] "India" "India"
## [3] "India" "Iran"
## [5] "India" "Yemen"
## [7] "India" "India"
## [9] "India" "India"
## [11] "Pakistan" "Russia"
## [13] "Russia" "Ukraine"
## [15] "Indonesia" "United States of America"
## [17] "Russia" "India"
## [19] "Pakistan" "Pakistan"
## [21] "Iran" "Indonesia"
## [23] "Russia" "India"
## [25] "India" "Pakistan"
## [27] "South Korea" "South Korea"
## [29] "South Korea" "South Korea"
## [31] "Indonesia" "India"
## [33] "Russia" "Taiwan"
## [35] "Taiwan" "Netherlands"
## [37] "Russia" "United Kingdom"
## [39] "India" "Indonesia"
## [41] "Tonga" "United Kingdom"
## [43] "Honduras" "Saudi Arabia"
## [45] "Canada" "Lithuania"
## [47] "India" "India"
## [49] "Turkey" "Turkey"
## [51] "United States of America" "United States of America"
## [53] "Australia" "India"
## [55] "India" "United States of America"
## [57] "Colombia" "United States of America"
## [59] "Australia" "United States of America"
## [61] "United States of America" "Mexico"
## [63] "Nigeria" "Bulgaria"
## [65] "Germany" "Venezuela"
## [67] "United States of America" "United Kingdom"
## [69] "United States of America" "Slovakia"
## [71] "United States of America" "United States of America"
## [73] "India" NA
## [75] NA "China"
## [77] "United States of America" "India"
## [79] "Romania" "Cuba"
## [81] "United States of America" "Philippines"
## [83] "United States of America" "Italy"
## [85] "Mexico" "United Kingdom"
## [87] NA "Mexico"
## [89] "United States of America" "Australia"
## [91] "United States of America" "Mexico"
## [93] "Ukraine" "Ukraine"
## [95] "Uruguay" "United Kingdom"
## [97] NA "United States of America"
## [99] "Mexico" "Denmark"
## [101] "United States of America" "United States of America"
## [103] "United States of America" "Guatemala"
## [105] "United States of America" "United States of America"
## [107] "The Bahamas" NA
## [109] "United States of America" "United States of America"
## [111] "United States of America" "Turkey"
## [113] "India" "Thailand"
## [115] "China" "Egypt"
## [117] "Ukraine" "United States of America"
## [119] "Nigeria" "United States of America"
## [121] "United States of America" "Nigeria"
## [123] "United States of America" "United Kingdom"
## [125] "United States of America" "India"
## [127] "India" "India"
## [129] "India" "United Republic of Tanzania"
## [131] "Nigeria" "Democratic Republic of the Congo"
## [133] "Ghana" "Chile"
## [135] "India" "United States of America"
## [137] "United Kingdom" "India"
## [139] "Turkey" "United States of America"
## [141] "United States of America" "United States of America"
## [143] "Mexico" "Mexico"
## [145] "United States of America" NA
## [147] "Syria" "United Kingdom"
## [149] NA "Jamaica"
## [151] "Myanmar" "United States of America"
## [153] "Ecuador" "United States of America"
## [155] "United States of America" "Nigeria"
## [157] "Democratic Republic of the Congo" "United States of America"
## [159] "United States of America" "Mexico"
## [161] "Italy" "Australia"
## [163] "Germany" "Ukraine"
## [165] "Australia" "Algeria"
## [167] "India" "Uruguay"
## [169] NA "Indonesia"
## [171] "South Africa" "Dominican Republic"
## [173] "Mexico" "India"
## [175] "United States of America" "Bangladesh"
## [177] "United States of America" "Australia"
## [179] "Russia" "Russia"
## [181] "Honduras" "United States of America"
## [183] "United States of America" "Denmark"
## [185] "United States of America" "United States of America"
## [187] "Russia" NA
## [189] "Mexico" "United States of America"
## [191] "Mexico" NA
## [193] NA "Indonesia"
## [195] "India" "Mexico"
## [197] "Montenegro" "United States of America"
## [199] "Italy" "China"
## [201] "Denmark" "Thailand"
## [203] "Zimbabwe" "India"
## [205] "India" "Cuba"
## [207] "Mexico" "Canada"
## [209] "Canada" "Canada"
## [211] "Canada" "Canada"
## [213] "China" "Canada"
## [215] "Canada" "Canada"
## [217] "Canada" "Canada"
## [219] "Canada" "Canada"
## [221] "Canada" NA
## [223] "Canada" "Canada"
## [225] "Canada" "Canada"
## [227] "Nigeria" "North Macedonia"
## [229] "India" "India"
## [231] "North Korea" "Poland"
## [233] "Kazakhstan" "China"
## [235] "China" "Russia"
## [237] "Iraq" "Russia"
## [239] "Russia" "Finland"
## [241] "France" "France"
## [243] "France" "France"
## [245] "France" "Japan"
## [247] "Japan" "Japan"
## [249] "Japan" "Norway"
## [251] "Brazil" "Brazil"
## [253] "Brazil" "Brazil"
## [255] "Brazil" "Portugal"
## [257] "Brazil" "Brazil"
## [259] "Brazil" "Brazil"
## [261] "Brazil" "Spain"
## [263] "Spain" "Spain"
## [265] "Sweden" "Sweden"
## [267] "Indonesia" "Argentina"
## [269] "Brazil" NA
## [271] "Egypt" "India"
## [273] "India" "Iran"
## [275] "Libya" "Saudi Arabia"
## [277] "Uruguay" "Kazakhstan"
## [279] "United States of America" "France"
## [281] "Netherlands" "Djibouti"
## [283] "Chile" "Turkey"
## [285] NA NA
## [287] "India" "Thailand"
## [289] "Fiji" "Mexico"
## [291] "Italy" "Jordan"
## [293] "Chile" "Mauritania"
## [295] "Colombia" "Argentina"
## [297] "Ivory Coast" "Belize"
## [299] "United States of America" "Netherlands"
## [301] "Peru" "United Republic of Tanzania"
## [303] "United States of America" "East Timor"
## [305] "Russia" "New Zealand"
## [307] NA "United States of America"
## [309] "Armenia" "United States of America"
## [311] NA "United States of America"
## [313] "Democratic Republic of the Congo" "Cameroon"
## [315] "United States of America" "United States of America"
## [317] "Comoros" "Australia"
## [319] "Solomon Islands" "Ukraine"
## [321] "Russia" "United States of America"
## [323] "Nauru" "Chile"
## [325] "Argentina" "China"
## [327] "China" "United Kingdom"
## [329] NA "Togo"
## [331] "Ukraine" NA
## [333] NA NA
## [335] "Mozambique" "Mauritius"
## [337] "Italy" "Mauritania"
## [339] "China" "Mauritania"
## [341] "France" NA
## [343] "United States of America" "Ireland"
## [345] "Indonesia" NA
## [347] "Pakistan" "United Kingdom"
## [349] "Brazil" "New Zealand"
## [351] NA "France"
## [353] "Italy" "Sudan"
## [355] "Cabo Verde" "Morocco"
## [357] "New Zealand" "Cambodia"
## [359] "Argentina" "Argentina"
## [361] "United States of America" NA
## [363] "Cabo Verde" "Bosnia and Herzegovina"
## [365] "Saint Kitts and Nevis" "Argentina"
## [367] "United States of America" "Bolivia"
## [369] NA "Yemen"
## [371] "China" NA
## [373] "Morocco" "Kiribati"
## [375] "Kazakhstan" "China"
## [377] "Australia" "Argentina"
## [379] "Tunisia" "United States of America"
## [381] "Mongolia" "China"
## [383] "Vanuatu" "Papua New Guinea"
## [385] "France" "China"
## [387] NA "Mexico"
## [389] "France" "Russia"
## [391] "China" "Azerbaijan"
## [393] "Cyprus" "Brazil"
## [395] "Finland" "Luxembourg"
## [397] "Italy" "South Korea"
## [399] "Brazil" "Brazil"
## [401] "Spain" "Argentina"
## [403] "Sweden" "Indonesia"
## [405] "Central African Republic" "China"
## [407] NA "Lithuania"
## [409] "New Zealand" "Paraguay"
## [411] "Peru" "Russia"
## [413] "Kosovo" "Belgium"
## [415] "Mexico" "China"
## [417] "India" "Turkey"
## [419] "China" "The Bahamas"
## [421] "Libya" "Canada"
## [423] "Qatar" "United States of America"
## [425] "Antigua and Barbuda" "Samoa"
## [427] "Egypt" NA
## [429] "United States of America" NA
## [431] "Gambia" "Burundi"
## [433] "Malawi" NA
## [435] "Iran" "Iraq"
## [437] "Mexico" "Sri Lanka"
## [439] NA "Thailand"
## [441] "United States of America" "Mexico"
## [443] "China" "United States of America"
## [445] "Vietnam" "United States of America"
## [447] "Syria" "Panama"
## [449] "Ethiopia" "Indonesia"
## [451] "United States of America" "Uganda"
## [453] "Democratic Republic of the Congo" NA
## [455] "India" "Poland"
## [457] NA "Russia"
## [459] "Mexico" "China"
## [461] "United States of America" "Mexico"
## [463] "United States of America" "United States of America"
## [465] "United States of America" "United States of America"
## [467] "Malaysia" "Rwanda"
## [469] "Poland" "Russia"
## [471] "China" "Philippines"
## [473] "United States of America" "United Kingdom"
## [475] "Slovenia" "India"
## [477] "Argentina" "Mexico"
## [479] "United Kingdom" "United States of America"
## [481] "Oman" "Tunisia"
## [483] "France" "Mexico"
## [485] "China" "Mexico"
## [487] "United States of America" "United States of America"
## [489] "United States of America" "United States of America"
## [491] "Suriname" "Russia"
## [493] "Malaysia" "Nigeria"
## [495] "Australia" "United States of America"
## [497] "United States of America" "United States of America"
## [499] "United Kingdom" "Liberia"
## [501] "India" "El Salvador"
## [503] "United States of America" "United States of America"
## [505] "United States of America" "Cuba"
## [507] "Mexico" "United States of America"
## [509] "Dominican Republic" "Russia"
## [511] "United States of America" "Iran"
## [513] "Mexico" "United States of America"
## [515] "United States of America" "Russia"
## [517] "Saint Lucia" "Namibia"
## [519] "Canada" "Canada"
## [521] "Canada" "Canada"
## [523] "Canada" "Canada"
## [525] "Canada" "Canada"
## [527] "Canada" NA
## [529] "China" "United States of America"
## [531] "United States of America" "China"
## [533] "China" "China"
## [535] "China" "Germany"
## [537] "Morocco" "France"
## [539] "Mexico" "France"
## [541] "Germany" "South Korea"
## [543] NA "Norway"
## [545] "Norway" "Brazil"
## [547] NA "Brazil"
## [549] "Brazil" "Spain"
## [551] "Spain" "Spain"
## [553] "Sweden" "Sweden"
## [555] "United Arab Emirates" "Algeria"
## [557] "Egypt" "Finland"
## [559] "Guyana" "India"
## [561] "Kenya" NA
## [563] "Mexico" "Pakistan"
## [565] "Belgium" "Nigeria"
## [567] "United States of America" "United States of America"
## [569] "Turkey" "United Kingdom"
## [571] "Italy" NA
## [573] "India" "Australia"
## [575] "United States of America" "Morocco"
## [577] "Germany" "Turkey"
## [579] "New Zealand" "India"
## [581] "United States of America" "Ukraine"
## [583] "Sudan" "United States of America"
## [585] "Belarus" NA
## [587] "United States of America" "United States of America"
## [589] "China" "United Arab Emirates"
## [591] "United States of America" "United States of America"
## [593] "Germany" "Costa Rica"
## [595] "China" "United States of America"
## [597] "Switzerland" "France"
## [599] NA "Portugal"
## [601] "Brazil" "Sweden"
## [603] "Indonesia" "Argentina"
## [605] "United Kingdom" "China"
## [607] "China" "Greece"
## [609] "Georgia" "Portugal"
## [611] "Taiwan" "Nigeria"
## [613] "Ecuador" "Pakistan"
## [615] "Moldova" "Peru"
## [617] "Canada" "Lithuania"
## [619] "China" "Mexico"
## [621] NA "Iran"
## [623] "Zimbabwe" "Kuwait"
## [625] "Canada" "Afghanistan"
## [627] NA "Mexico"
## [629] "Ghana" "Ethiopia"
## [631] "Yemen" "Australia"
## [633] "Kazakhstan" "Algeria"
## [635] "Syria" "India"
## [637] "Chile" "Turkmenistan"
## [639] "Eritrea" "Paraguay"
## [641] "United Kingdom" "Republic of Serbia"
## [643] "Lebanon" "India"
## [645] "Mali" "United States of America"
## [647] "Australia" "United States of America"
## [649] "United States of America" "United States of America"
## [651] "Zimbabwe" "Brunei"
## [653] "China" "Chile"
## [655] "India" NA
## [657] "New Zealand" "Guinea"
## [659] "United States of America" "Colombia"
## [661] "Benin" "Argentina"
## [663] "Colombia" "Mexico"
## [665] "Mexico" "Bangladesh"
## [667] NA "South Africa"
## [669] "Democratic Republic of the Congo" "Morocco"
## [671] "Democratic Republic of the Congo" NA
## [673] "North Korea" "Kyrgyzstan"
## [675] "Botswana" "Mexico"
## [677] "United Kingdom" "Guatemala"
## [679] NA "Vietnam"
## [681] "Cuba" "Egypt"
## [683] "Saudi Arabia" "Nigeria"
## [685] "China" NA
## [687] "Nepal" "Angola"
## [689] "Russia" "Pakistan"
## [691] "Malawi" "Nigeria"
## [693] "Bolivia" "Zambia"
## [695] "Egypt" "India"
## [697] "Venezuela" "Colombia"
## [699] "United States of America" "Nicaragua"
## [701] "Iran" "United States of America"
## [703] "Mexico" "Malta"
## [705] "Kenya" "Lesotho"
## [707] "United States of America" "Pakistan"
## [709] "Uruguay" "Mexico"
## [711] "The Bahamas" "Chad"
## [713] "Niger" NA
## [715] "Angola" "United States of America"
## [717] NA "Burkina Faso"
## [719] "Haiti" "Mexico"
## [721] "United States of America" "Australia"
## [723] "South Africa" "Chile"
## [725] "Cambodia" "India"
## [727] "Papua New Guinea" "Panama"
## [729] "Chile" "United States of America"
## [731] "Myanmar" "Latvia"
## [733] "Yemen" "Iraq"
## [735] NA "Vietnam"
## [737] "Greece" "Bulgaria"
## [739] NA "Fiji"
## [741] "Iran" "Mexico"
## [743] "Honduras" "Iran"
## [745] "Albania" "United States of America"
## [747] "Mexico" "Estonia"
## [749] "Israel" NA
## [751] "Madagascar" "United States of America"
## [753] "Venezuela" "Russia"
## [755] "Laos" "Bolivia"
## [757] NA NA
## [759] "Canada" "Canada"
## [761] "Canada" "Canada"
## [763] "Canada" "Canada"
## [765] "Canada" "Canada"
## [767] "Canada" "Canada"
## [769] "Canada" "Croatia"
## [771] "United Republic of Tanzania" "Iceland"
## [773] "Russia" "Russia"
## [775] "France" "Canada"
## [777] NA "Japan"
## [779] NA "Japan"
## [781] "Norway" NA
## [783] NA "Brazil"
## [785] "Brazil" "Brazil"
## [787] "Brazil" "Brazil"
## [789] "Brazil" "Brazil"
## [791] "Brazil" "Argentina"
## [793] "Brazil" "Cameroon"
## [795] NA "Jordan"
## [797] NA "Mexico"
## [799] "Guinea-Bissau" "Philippines"
## [801] "Seychelles" "Senegal"
## [803] "Sudan" "Uzbekistan"
## [805] "Denmark" "Romania"
## [807] "Hungary" "China"
## [809] "United States of America" "United States of America"
## [811] "Ireland" "United States of America"
## [813] "Germany" "Malaysia"
## [815] "United States of America" "United States of America"
## [817] "United States of America" "Germany"
## [819] NA "United States of America"
## [821] "United States of America" "United States of America"
## [823] "Poland" "Switzerland"
## [825] "Belgium" "Germany"
## [827] "Spain" "Czechia"
## [829] NA "Taiwan"
## [831] "Netherlands" "Singapore"
## [833] "United Kingdom" NA
## [835] "United States of America" "United States of America"
## [837] "China" "Colombia"
## [839] "India" "United States of America"
## [841] "United States of America" "Egypt"
## [843] "Morocco" "Venezuela"
## [845] "South Africa" "China"
## [847] "India" "United States of America"
## [849] "United States of America" "Thailand"
## [851] "United Arab Emirates" "United States of America"
## [853] "Argentina" "United States of America"
## [855] "United States of America" NA
## [857] "South Africa" NA
## [859] "United States of America" "Italy"
## [861] "Australia" "Mexico"
## [863] "Philippines" "Kenya"
## [865] NA "United States of America"
## [867] "Saudi Arabia" "Chile"
## [869] "United States of America" "United States of America"
## [871] "China" "Russia"
## [873] "Canada" "Australia"
## [875] "Finland" "France"
## [877] "Germany" "Austria"
## [879] "Germany" "Italy"
## [881] "Japan" "South Korea"
## [883] "Norway" "Brazil"
## [885] "Brazil" "Brazil"
## [887] "Spain" "United States of America"
## [889] "Sweden" "Indonesia"
## [891] "Greece" NA
## [893] "Thailand"#check how many airports per country
airport_count_by_country <- airports_with_country %>%
group_by(sovereignt) %>%
summarise(airport_count = n())
print(airport_count_by_country)## Simple feature collection with 171 features and 2 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -175.1356 ymin: -53.78147 xmax: 179.1954 ymax: 78.24672
## Geodetic CRS: WGS 84
## # A tibble: 171 × 3
## sovereignt airport_count geometry
## <chr> <int> <GEOMETRY [°]>
## 1 Afghanistan 1 POINT (69.21007 34.5634)
## 2 Albania 1 POINT (19.71503 41.42085)
## 3 Algeria 3 MULTIPOINT ((-0.606797 35.62027), (3.21207…
## 4 Angola 2 MULTIPOINT ((13.2348 -8.848313), (15.74976…
## 5 Antigua and Barbuda 1 POINT (-61.79237 17.14036)
## 6 Argentina 13 MULTIPOINT ((-65.09377 -24.3861), (-65.478…
## 7 Armenia 1 POINT (44.40006 40.15237)
## 8 Australia 15 MULTIPOINT ((115.9742 -31.94113), (133.902…
## 9 Austria 1 POINT (16.56077 48.11976)
## 10 Azerbaijan 1 POINT (50.04984 40.46275)
## # ℹ 161 more rowsairports_cleaned <- airports_with_country %>%
filter(!is.na(sovereignt))# First, let's see what columns we have in airports_cleaned
print("Columns in airports_cleaned:")## [1] "Columns in airports_cleaned:"names(airports_cleaned)## [1] "scalerank.x" "featurecla.x" "type.x" "name.x" "abbrev.x"
## [6] "location" "gps_code" "iata_code" "wikipedia" "natlscale"
## [11] "comments" "wikidataid.x" "name_ar.x" "name_bn.x" "name_de.x"
## [16] "name_en.x" "name_es.x" "name_fr.x" "name_el.x" "name_hi.x"
## [21] "name_hu.x" "name_id.x" "name_it.x" "name_ja.x" "name_ko.x"
## [26] "name_nl.x" "name_pl.x" "name_pt.x" "name_ru.x" "name_sv.x"
## [31] "name_tr.x" "name_vi.x" "name_zh.x" "wdid_score" "ne_id.x"
## [36] "name_fa.x" "name_he.x" "name_uk.x" "name_ur.x" "name_zht.x"
## [41] "featurecla.y" "scalerank.y" "labelrank" "sovereignt" "sov_a3"
## [46] "adm0_dif" "level" "type.y" "tlc" "admin"
## [51] "adm0_a3" "geou_dif" "geounit" "gu_a3" "su_dif"
## [56] "subunit" "su_a3" "brk_diff" "name.y" "name_long"
## [61] "brk_a3" "brk_name" "brk_group" "abbrev.y" "postal"
## [66] "formal_en" "formal_fr" "name_ciawf" "note_adm0" "note_brk"
## [71] "name_sort" "name_alt" "mapcolor7" "mapcolor8" "mapcolor9"
## [76] "mapcolor13" "pop_est" "pop_rank" "pop_year" "gdp_md"
## [81] "gdp_year" "economy" "income_grp" "fips_10" "iso_a2"
## [86] "iso_a2_eh" "iso_a3" "iso_a3_eh" "iso_n3" "iso_n3_eh"
## [91] "un_a3" "wb_a2" "wb_a3" "woe_id" "woe_id_eh"
## [96] "woe_note" "adm0_iso" "adm0_diff" "adm0_tlc" "adm0_a3_us"
## [101] "adm0_a3_fr" "adm0_a3_ru" "adm0_a3_es" "adm0_a3_cn" "adm0_a3_tw"
## [106] "adm0_a3_in" "adm0_a3_np" "adm0_a3_pk" "adm0_a3_de" "adm0_a3_gb"
## [111] "adm0_a3_br" "adm0_a3_il" "adm0_a3_ps" "adm0_a3_sa" "adm0_a3_eg"
## [116] "adm0_a3_ma" "adm0_a3_pt" "adm0_a3_ar" "adm0_a3_jp" "adm0_a3_ko"
## [121] "adm0_a3_vn" "adm0_a3_tr" "adm0_a3_id" "adm0_a3_pl" "adm0_a3_gr"
## [126] "adm0_a3_it" "adm0_a3_nl" "adm0_a3_se" "adm0_a3_bd" "adm0_a3_ua"
## [131] "adm0_a3_un" "adm0_a3_wb" "continent" "region_un" "subregion"
## [136] "region_wb" "name_len" "long_len" "abbrev_len" "tiny"
## [141] "homepart" "min_zoom" "min_label" "max_label" "label_x"
## [146] "label_y" "ne_id.y" "wikidataid.y" "name_ar.y" "name_bn.y"
## [151] "name_de.y" "name_en.y" "name_es.y" "name_fa.y" "name_fr.y"
## [156] "name_el.y" "name_he.y" "name_hi.y" "name_hu.y" "name_id.y"
## [161] "name_it.y" "name_ja.y" "name_ko.y" "name_nl.y" "name_pl.y"
## [166] "name_pt.y" "name_ru.y" "name_sv.y" "name_tr.y" "name_uk.y"
## [171] "name_ur.y" "name_vi.y" "name_zh.y" "name_zht.y" "fclass_iso"
## [176] "tlc_diff" "fclass_tlc" "fclass_us" "fclass_fr" "fclass_ru"
## [181] "fclass_es" "fclass_cn" "fclass_tw" "fclass_in" "fclass_np"
## [186] "fclass_pk" "fclass_de" "fclass_gb" "fclass_br" "fclass_il"
## [191] "fclass_ps" "fclass_sa" "fclass_eg" "fclass_ma" "fclass_pt"
## [196] "fclass_ar" "fclass_jp" "fclass_ko" "fclass_vn" "fclass_tr"
## [201] "fclass_id" "fclass_pl" "fclass_gr" "fclass_it" "fclass_nl"
## [206] "fclass_se" "fclass_bd" "fclass_ua" "geometry"# Now let's check step by step
pop_airport_pairs <- top20_pop %>%
group_by(sovereignt) %>%
mutate(matching_airports = list(airports_cleaned[airports_cleaned$sovereignt == first(sovereignt), ])) %>%
unnest(matching_airports, names_repair = "unique_quiet", names_sep = "_airport")
# Check what columns we actually have
print("Columns in pop_airport_pairs:")## [1] "Columns in pop_airport_pairs:"names(pop_airport_pairs)## [1] "SCALERANK"
## [2] "NATSCALE"
## [3] "LABELRANK"
## [4] "FEATURECLA"
## [5] "NAME"
## [6] "NAMEPAR"
## [7] "NAMEALT"
## [8] "NAMEASCII"
## [9] "ADM0CAP"
## [10] "CAPIN"
## [11] "WORLDCITY"
## [12] "MEGACITY"
## [13] "SOV0NAME"
## [14] "SOV_A3"
## [15] "ADM0NAME"
## [16] "ADM0_A3"
## [17] "ADM1NAME"
## [18] "ISO_A2"
## [19] "NOTE"
## [20] "LATITUDE"
## [21] "LONGITUDE"
## [22] "POP_MAX"
## [23] "POP_MIN"
## [24] "POP_OTHER"
## [25] "RANK_MAX"
## [26] "RANK_MIN"
## [27] "MEGANAME"
## [28] "LS_NAME"
## [29] "MAX_POP10"
## [30] "MAX_POP20"
## [31] "MAX_POP50"
## [32] "MAX_POP300"
## [33] "MAX_POP310"
## [34] "MAX_NATSCA"
## [35] "MIN_AREAKM"
## [36] "MAX_AREAKM"
## [37] "MIN_AREAMI"
## [38] "MAX_AREAMI"
## [39] "MIN_PERKM"
## [40] "MAX_PERKM"
## [41] "MIN_PERMI"
## [42] "MAX_PERMI"
## [43] "MIN_BBXMIN"
## [44] "MAX_BBXMIN"
## [45] "MIN_BBXMAX"
## [46] "MAX_BBXMAX"
## [47] "MIN_BBYMIN"
## [48] "MAX_BBYMIN"
## [49] "MIN_BBYMAX"
## [50] "MAX_BBYMAX"
## [51] "MEAN_BBXC"
## [52] "MEAN_BBYC"
## [53] "TIMEZONE"
## [54] "UN_FID"
## [55] "POP1950"
## [56] "POP1955"
## [57] "POP1960"
## [58] "POP1965"
## [59] "POP1970"
## [60] "POP1975"
## [61] "POP1980"
## [62] "POP1985"
## [63] "POP1990"
## [64] "POP1995"
## [65] "POP2000"
## [66] "POP2005"
## [67] "POP2010"
## [68] "POP2015"
## [69] "POP2020"
## [70] "POP2025"
## [71] "POP2050"
## [72] "MIN_ZOOM"
## [73] "WIKIDATAID"
## [74] "WOF_ID"
## [75] "CAPALT"
## [76] "NAME_EN"
## [77] "NAME_DE"
## [78] "NAME_ES"
## [79] "NAME_FR"
## [80] "NAME_PT"
## [81] "NAME_RU"
## [82] "NAME_ZH"
## [83] "LABEL"
## [84] "NAME_AR"
## [85] "NAME_BN"
## [86] "NAME_EL"
## [87] "NAME_HI"
## [88] "NAME_HU"
## [89] "NAME_ID"
## [90] "NAME_IT"
## [91] "NAME_JA"
## [92] "NAME_KO"
## [93] "NAME_NL"
## [94] "NAME_PL"
## [95] "NAME_SV"
## [96] "NAME_TR"
## [97] "NAME_VI"
## [98] "NE_ID"
## [99] "NAME_FA"
## [100] "NAME_HE"
## [101] "NAME_UK"
## [102] "NAME_UR"
## [103] "NAME_ZHT"
## [104] "GEONAMESID"
## [105] "FCLASS_ISO"
## [106] "FCLASS_US"
## [107] "FCLASS_FR"
## [108] "FCLASS_RU"
## [109] "FCLASS_ES"
## [110] "FCLASS_CN"
## [111] "FCLASS_TW"
## [112] "FCLASS_IN"
## [113] "FCLASS_NP"
## [114] "FCLASS_PK"
## [115] "FCLASS_DE"
## [116] "FCLASS_GB"
## [117] "FCLASS_BR"
## [118] "FCLASS_IL"
## [119] "FCLASS_PS"
## [120] "FCLASS_SA"
## [121] "FCLASS_EG"
## [122] "FCLASS_MA"
## [123] "FCLASS_PT"
## [124] "FCLASS_AR"
## [125] "FCLASS_JP"
## [126] "FCLASS_KO"
## [127] "FCLASS_VN"
## [128] "FCLASS_TR"
## [129] "FCLASS_ID"
## [130] "FCLASS_PL"
## [131] "FCLASS_GR"
## [132] "FCLASS_IT"
## [133] "FCLASS_NL"
## [134] "FCLASS_SE"
## [135] "FCLASS_BD"
## [136] "FCLASS_UA"
## [137] "FCLASS_TLC"
## [138] "sovereignt"
## [139] "continent"
## [140] "geometry"
## [141] "matching_airports_airportscalerank.x"
## [142] "matching_airports_airportfeaturecla.x"
## [143] "matching_airports_airporttype.x"
## [144] "matching_airports_airportname.x"
## [145] "matching_airports_airportabbrev.x"
## [146] "matching_airports_airportlocation"
## [147] "matching_airports_airportgps_code"
## [148] "matching_airports_airportiata_code"
## [149] "matching_airports_airportwikipedia"
## [150] "matching_airports_airportnatlscale"
## [151] "matching_airports_airportcomments"
## [152] "matching_airports_airportwikidataid.x"
## [153] "matching_airports_airportname_ar.x"
## [154] "matching_airports_airportname_bn.x"
## [155] "matching_airports_airportname_de.x"
## [156] "matching_airports_airportname_en.x"
## [157] "matching_airports_airportname_es.x"
## [158] "matching_airports_airportname_fr.x"
## [159] "matching_airports_airportname_el.x"
## [160] "matching_airports_airportname_hi.x"
## [161] "matching_airports_airportname_hu.x"
## [162] "matching_airports_airportname_id.x"
## [163] "matching_airports_airportname_it.x"
## [164] "matching_airports_airportname_ja.x"
## [165] "matching_airports_airportname_ko.x"
## [166] "matching_airports_airportname_nl.x"
## [167] "matching_airports_airportname_pl.x"
## [168] "matching_airports_airportname_pt.x"
## [169] "matching_airports_airportname_ru.x"
## [170] "matching_airports_airportname_sv.x"
## [171] "matching_airports_airportname_tr.x"
## [172] "matching_airports_airportname_vi.x"
## [173] "matching_airports_airportname_zh.x"
## [174] "matching_airports_airportwdid_score"
## [175] "matching_airports_airportne_id.x"
## [176] "matching_airports_airportname_fa.x"
## [177] "matching_airports_airportname_he.x"
## [178] "matching_airports_airportname_uk.x"
## [179] "matching_airports_airportname_ur.x"
## [180] "matching_airports_airportname_zht.x"
## [181] "matching_airports_airportfeaturecla.y"
## [182] "matching_airports_airportscalerank.y"
## [183] "matching_airports_airportlabelrank"
## [184] "matching_airports_airportsovereignt"
## [185] "matching_airports_airportsov_a3"
## [186] "matching_airports_airportadm0_dif"
## [187] "matching_airports_airportlevel"
## [188] "matching_airports_airporttype.y"
## [189] "matching_airports_airporttlc"
## [190] "matching_airports_airportadmin"
## [191] "matching_airports_airportadm0_a3"
## [192] "matching_airports_airportgeou_dif"
## [193] "matching_airports_airportgeounit"
## [194] "matching_airports_airportgu_a3"
## [195] "matching_airports_airportsu_dif"
## [196] "matching_airports_airportsubunit"
## [197] "matching_airports_airportsu_a3"
## [198] "matching_airports_airportbrk_diff"
## [199] "matching_airports_airportname.y"
## [200] "matching_airports_airportname_long"
## [201] "matching_airports_airportbrk_a3"
## [202] "matching_airports_airportbrk_name"
## [203] "matching_airports_airportbrk_group"
## [204] "matching_airports_airportabbrev.y"
## [205] "matching_airports_airportpostal"
## [206] "matching_airports_airportformal_en"
## [207] "matching_airports_airportformal_fr"
## [208] "matching_airports_airportname_ciawf"
## [209] "matching_airports_airportnote_adm0"
## [210] "matching_airports_airportnote_brk"
## [211] "matching_airports_airportname_sort"
## [212] "matching_airports_airportname_alt"
## [213] "matching_airports_airportmapcolor7"
## [214] "matching_airports_airportmapcolor8"
## [215] "matching_airports_airportmapcolor9"
## [216] "matching_airports_airportmapcolor13"
## [217] "matching_airports_airportpop_est"
## [218] "matching_airports_airportpop_rank"
## [219] "matching_airports_airportpop_year"
## [220] "matching_airports_airportgdp_md"
## [221] "matching_airports_airportgdp_year"
## [222] "matching_airports_airporteconomy"
## [223] "matching_airports_airportincome_grp"
## [224] "matching_airports_airportfips_10"
## [225] "matching_airports_airportiso_a2"
## [226] "matching_airports_airportiso_a2_eh"
## [227] "matching_airports_airportiso_a3"
## [228] "matching_airports_airportiso_a3_eh"
## [229] "matching_airports_airportiso_n3"
## [230] "matching_airports_airportiso_n3_eh"
## [231] "matching_airports_airportun_a3"
## [232] "matching_airports_airportwb_a2"
## [233] "matching_airports_airportwb_a3"
## [234] "matching_airports_airportwoe_id"
## [235] "matching_airports_airportwoe_id_eh"
## [236] "matching_airports_airportwoe_note"
## [237] "matching_airports_airportadm0_iso"
## [238] "matching_airports_airportadm0_diff"
## [239] "matching_airports_airportadm0_tlc"
## [240] "matching_airports_airportadm0_a3_us"
## [241] "matching_airports_airportadm0_a3_fr"
## [242] "matching_airports_airportadm0_a3_ru"
## [243] "matching_airports_airportadm0_a3_es"
## [244] "matching_airports_airportadm0_a3_cn"
## [245] "matching_airports_airportadm0_a3_tw"
## [246] "matching_airports_airportadm0_a3_in"
## [247] "matching_airports_airportadm0_a3_np"
## [248] "matching_airports_airportadm0_a3_pk"
## [249] "matching_airports_airportadm0_a3_de"
## [250] "matching_airports_airportadm0_a3_gb"
## [251] "matching_airports_airportadm0_a3_br"
## [252] "matching_airports_airportadm0_a3_il"
## [253] "matching_airports_airportadm0_a3_ps"
## [254] "matching_airports_airportadm0_a3_sa"
## [255] "matching_airports_airportadm0_a3_eg"
## [256] "matching_airports_airportadm0_a3_ma"
## [257] "matching_airports_airportadm0_a3_pt"
## [258] "matching_airports_airportadm0_a3_ar"
## [259] "matching_airports_airportadm0_a3_jp"
## [260] "matching_airports_airportadm0_a3_ko"
## [261] "matching_airports_airportadm0_a3_vn"
## [262] "matching_airports_airportadm0_a3_tr"
## [263] "matching_airports_airportadm0_a3_id"
## [264] "matching_airports_airportadm0_a3_pl"
## [265] "matching_airports_airportadm0_a3_gr"
## [266] "matching_airports_airportadm0_a3_it"
## [267] "matching_airports_airportadm0_a3_nl"
## [268] "matching_airports_airportadm0_a3_se"
## [269] "matching_airports_airportadm0_a3_bd"
## [270] "matching_airports_airportadm0_a3_ua"
## [271] "matching_airports_airportadm0_a3_un"
## [272] "matching_airports_airportadm0_a3_wb"
## [273] "matching_airports_airportcontinent"
## [274] "matching_airports_airportregion_un"
## [275] "matching_airports_airportsubregion"
## [276] "matching_airports_airportregion_wb"
## [277] "matching_airports_airportname_len"
## [278] "matching_airports_airportlong_len"
## [279] "matching_airports_airportabbrev_len"
## [280] "matching_airports_airporttiny"
## [281] "matching_airports_airporthomepart"
## [282] "matching_airports_airportmin_zoom"
## [283] "matching_airports_airportmin_label"
## [284] "matching_airports_airportmax_label"
## [285] "matching_airports_airportlabel_x"
## [286] "matching_airports_airportlabel_y"
## [287] "matching_airports_airportne_id.y"
## [288] "matching_airports_airportwikidataid.y"
## [289] "matching_airports_airportname_ar.y"
## [290] "matching_airports_airportname_bn.y"
## [291] "matching_airports_airportname_de.y"
## [292] "matching_airports_airportname_en.y"
## [293] "matching_airports_airportname_es.y"
## [294] "matching_airports_airportname_fa.y"
## [295] "matching_airports_airportname_fr.y"
## [296] "matching_airports_airportname_el.y"
## [297] "matching_airports_airportname_he.y"
## [298] "matching_airports_airportname_hi.y"
## [299] "matching_airports_airportname_hu.y"
## [300] "matching_airports_airportname_id.y"
## [301] "matching_airports_airportname_it.y"
## [302] "matching_airports_airportname_ja.y"
## [303] "matching_airports_airportname_ko.y"
## [304] "matching_airports_airportname_nl.y"
## [305] "matching_airports_airportname_pl.y"
## [306] "matching_airports_airportname_pt.y"
## [307] "matching_airports_airportname_ru.y"
## [308] "matching_airports_airportname_sv.y"
## [309] "matching_airports_airportname_tr.y"
## [310] "matching_airports_airportname_uk.y"
## [311] "matching_airports_airportname_ur.y"
## [312] "matching_airports_airportname_vi.y"
## [313] "matching_airports_airportname_zh.y"
## [314] "matching_airports_airportname_zht.y"
## [315] "matching_airports_airportfclass_iso"
## [316] "matching_airports_airporttlc_diff"
## [317] "matching_airports_airportfclass_tlc"
## [318] "matching_airports_airportfclass_us"
## [319] "matching_airports_airportfclass_fr"
## [320] "matching_airports_airportfclass_ru"
## [321] "matching_airports_airportfclass_es"
## [322] "matching_airports_airportfclass_cn"
## [323] "matching_airports_airportfclass_tw"
## [324] "matching_airports_airportfclass_in"
## [325] "matching_airports_airportfclass_np"
## [326] "matching_airports_airportfclass_pk"
## [327] "matching_airports_airportfclass_de"
## [328] "matching_airports_airportfclass_gb"
## [329] "matching_airports_airportfclass_br"
## [330] "matching_airports_airportfclass_il"
## [331] "matching_airports_airportfclass_ps"
## [332] "matching_airports_airportfclass_sa"
## [333] "matching_airports_airportfclass_eg"
## [334] "matching_airports_airportfclass_ma"
## [335] "matching_airports_airportfclass_pt"
## [336] "matching_airports_airportfclass_ar"
## [337] "matching_airports_airportfclass_jp"
## [338] "matching_airports_airportfclass_ko"
## [339] "matching_airports_airportfclass_vn"
## [340] "matching_airports_airportfclass_tr"
## [341] "matching_airports_airportfclass_id"
## [342] "matching_airports_airportfclass_pl"
## [343] "matching_airports_airportfclass_gr"
## [344] "matching_airports_airportfclass_it"
## [345] "matching_airports_airportfclass_nl"
## [346] "matching_airports_airportfclass_se"
## [347] "matching_airports_airportfclass_bd"
## [348] "matching_airports_airportfclass_ua"
## [349] "matching_airports_airportgeometry"within_country_distances <- pop_airport_pairs %>%
select(
sovereignt,
city = NAME,
airport = matching_airports_airportname.x,
geometry,
geometry_airport = matching_airports_airportgeometry
) %>%
rowwise() %>% # Process each row individually
mutate(
distance_km = as.numeric(st_distance(geometry, geometry_airport))/1000
) %>%
ungroup() %>% # Remove rowwise grouping
arrange(sovereignt, distance_km)country_avg_distances <- within_country_distances %>%
group_by(sovereignt) %>%
summarise(
avg_distance_km = mean(distance_km)
) %>%
arrange(avg_distance_km)
# Let's look at the results
print(country_avg_distances)## Simple feature collection with 157 features and 2 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -175.2206 ymin: -41.29207 xmax: 174.7772 ymax: 64.14346
## Geodetic CRS: WGS 84
## # A tibble: 157 × 3
## sovereignt avg_distance_km geometry
## <chr> <dbl> <POINT [°]>
## 1 Saint Kitts and Nevis 1.05 (-62.71701 17.30203)
## 2 Iceland 1.38 (-21.93655 64.14346)
## 3 Chad 2.23 (15.0472 12.11504)
## 4 Burkina Faso 2.50 (-1.52667 12.37226)
## 5 Qatar 3.28 (51.53297 25.28656)
## 6 Mozambique 3.48 (32.58722 -25.95333)
## 7 Laos 3.50 (102.6 17.96669)
## 8 Senegal 3.51 (-17.47508 14.71778)
## 9 Vanuatu 3.53 (168.3166 -17.73335)
## 10 Turkmenistan 4.23 (58.3833 37.94999)
## # ℹ 147 more rows# Create continent-wise data
continent_distances <- country_avg_distances %>%
left_join(
pop_airport_pairs %>%
st_drop_geometry() %>% # Remove the geometry to make it a regular dataframe
select(sovereignt, continent) %>%
distinct(),
by = "sovereignt"
)
print(continent_distances)## Simple feature collection with 157 features and 3 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -175.2206 ymin: -41.29207 xmax: 174.7772 ymax: 64.14346
## Geodetic CRS: WGS 84
## # A tibble: 157 × 4
## sovereignt avg_distance_km geometry continent
## <chr> <dbl> <POINT [°]> <chr>
## 1 Saint Kitts and Nevis 1.05 (-62.71701 17.30203) North America
## 2 Iceland 1.38 (-21.93655 64.14346) Europe
## 3 Chad 2.23 (15.0472 12.11504) Africa
## 4 Burkina Faso 2.50 (-1.52667 12.37226) Africa
## 5 Qatar 3.28 (51.53297 25.28656) Asia
## 6 Mozambique 3.48 (32.58722 -25.95333) Africa
## 7 Laos 3.50 (102.6 17.96669) Asia
## 8 Senegal 3.51 (-17.47508 14.71778) Africa
## 9 Vanuatu 3.53 (168.3166 -17.73335) Oceania
## 10 Turkmenistan 4.23 (58.3833 37.94999) Asia
## # ℹ 147 more rowsggplot(continent_distances, aes(x = avg_distance_km)) +
geom_histogram(bins = 10, fill = "red", color = "black") +
facet_wrap(~continent) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 0, hjust = 0.5, size = 8),
axis.text.y = element_text(size = 8),
strip.text = element_text(size = 10), # Continent labels size
panel.grid.minor = element_blank(),
plot.title = element_text(size = 14) # Title size
) +
scale_x_continuous(
breaks = seq(0, 4000, by = 1000),
labels = seq(0, 4000, by = 1000)
) +
scale_y_continuous(
breaks = seq(0, 25, by = 5) # Set specific y-axis breaks
) +
labs(
x = "Average Distance (km)",
y = "Count",
title = "Distribution of Average Airport Distances by Continent"
)
## THE OUTLIERS IN THE HISTOGRAMS ARE PRESENT SINCE THERE ARE COUNTED THE TERRITORIES WHICH BELONG TO EUROPEAN COUNTRIES, FOR EXAMPLE QUEEN BEATRIZ INTERNATIONAL AIRPORT SO OUR FUNCTION CALCULATES DISTANCE BETWEEN AMSTERDAM AND DUTCH CARIBBEAN ISLAND OF ARUBA6 Prepare data
# Download and read the market coordinates data from a GitHub raw file
url <- "https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/Data/Part%202/MktCoords.xlsx"
temp_file <- tempfile(fileext = ".xlsx")
download.file(url, temp_file, mode = "wb")
market_data <- read_excel(temp_file)# Convert the market coordinates dataframe into a spatial sf object using longitude and latitude
market_sf <- st_as_sf(market_data, coords = c("longitude", "latitude"), crs = 4326)
print(market_sf)## Simple feature collection with 230 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -17.36603 ymin: -29.3173 xmax: 49.18333 ymax: 19.61667
## Geodetic CRS: WGS 84
## # A tibble: 230 × 4
## ctrycode mktcode market geometry
## * <chr> <dbl> <chr> <POINT [°]>
## 1 AO 200 Luanda (13.23444 -8.838333)
## 2 BJ 401 Cotonou (2.433333 6.366667)
## 3 BJ 402 Malanville (3.393031 11.86195)
## 4 BJ 403 Natitingou (1.388616 10.31291)
## 5 BJ 404 Parakou (2.616667 9.35)
## 6 BW 110 Gaborone (25.90874 -24.65411)
## 7 BF 405 Bobo Dioulasso (-4.291773 11.1781)
## 8 BF 406 Dedougou (-3.463617 12.45462)
## 9 BF 407 Fada Ngourma (0.360972 12.05481)
## 10 BF 408 Ouagadougou (-1.533864 12.36464)
## # ℹ 220 more rowsworld_map <- map_data("world")# Filter for Africa’s
africa_map <- subset(world_map, long > -25 & long < 60 & lat > -40 & lat < 40)
# Create plot
ggplot() +
geom_polygon(data = africa_map, aes(x = long, y = lat, group = group), fill = "gray90", color = "white") +
geom_sf(data = market_sf, color = "red", size = 2, alpha = 0.7) +
theme_minimal() +
labs(title = "Market Locations") + # Remove x and y labels
theme(
axis.title = element_blank(), # Remove axis titles
axis.text = element_blank(), # Remove axis text
axis.ticks = element_blank(), # Remove axis ticks
plot.title = element_text(hjust = 0.5, face = "bold")
)
# Download and read the price data from a GitHub raw file
price_file <- "https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/Data/Part%202/PriceMaster4GAMS.xlsx"
temp_file <- tempfile(fileext = ".xlsx")
download.file(price_file, temp_file, mode = "wb")
price_data <- read_excel(temp_file)
head(price_data)## # A tibble: 6 × 148
## mktcode country market crop `1` `2` `3` `4` `5` `6` `7`
## <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 200 Angola Luanda Maize NA NA NA NA NA NA NA
## 2 401 Benin Cotonou Maize 0.250 0.250 0.263 0.277 0.299 0.300 0.311
## 3 401 Benin Cotonou Rice 0.480 0.480 0.472 0.467 0.481 0.504 0.532
## 4 402 Benin Malanv… Maize 0.169 0.176 0.189 0.203 0.216 0.205 0.266
## 5 402 Benin Malanv… Sorg… 0.182 0.189 0.209 0.216 0.258 0.300 0.319
## 6 403 Benin Natiti… Maize 0.189 0.216 0.216 0.243 0.258 0.270 0.311
## # ℹ 137 more variables: `8` <dbl>, `9` <dbl>, `10` <dbl>, `11` <dbl>,
## # `12` <dbl>, `13` <dbl>, `14` <dbl>, `15` <dbl>, `16` <dbl>, `17` <dbl>,
## # `18` <dbl>, `19` <dbl>, `20` <dbl>, `21` <dbl>, `22` <dbl>, `23` <dbl>,
## # `24` <dbl>, `25` <dbl>, `26` <dbl>, `27` <dbl>, `28` <dbl>, `29` <dbl>,
## # `30` <dbl>, `31` <dbl>, `32` <dbl>, `33` <dbl>, `34` <dbl>, `35` <dbl>,
## # `36` <dbl>, `37` <dbl>, `38` <dbl>, `39` <dbl>, `40` <dbl>, `41` <dbl>,
## # `42` <dbl>, `43` <dbl>, `44` <dbl>, `45` <dbl>, `46` <dbl>, `47` <dbl>, …# Convert price columns to numeric format, ignoring the first 4 columns
price_data <- price_data %>%
mutate(across(5:ncol(.), as.numeric, .names = "num_{.col}"))To calculate average market prices, we have taken the average price across all columns, assuming they represent different years, months, or days. Additionally, we compute the average price per market code (mktcode), without differentiating by crop type. This is because, for our analysis, we are interested in understanding market-level price trends rather than specific crop variations.
# Compute the average price per market code by averaging across all price columns,
# assuming they represent different time periods (e.g., years, months, or days),
# without differentiating by crop type since it's not relevant for our analysis.
average_price <- price_data %>%
rowwise() %>%
mutate(avg_price = mean(c_across(5:ncol(.)), na.rm = TRUE)) %>%
ungroup() %>%
select(mktcode, avg_price) %>%
group_by(mktcode) %>%
summarise(avg_price = mean(avg_price, na.rm = TRUE))
# Print first few rows to verify the computed average prices
print(head(average_price))## # A tibble: 6 × 2
## mktcode avg_price
## <dbl> <dbl>
## 1 110 0.262
## 2 121 0.838
## 3 122 0.897
## 4 123 0.986
## 5 124 0.912
## 6 130 0.305# Merge market coordinates with the computed average prices
# using "mktcode" as the common key to link both datasets.
merged_data <- market_sf %>%
left_join(average_price, by = "mktcode")
# Print the first few rows to verify that the merge was successful
print(head(merged_data))## Simple feature collection with 6 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1.388616 ymin: -24.65411 xmax: 25.90874 ymax: 11.86195
## Geodetic CRS: WGS 84
## # A tibble: 6 × 5
## ctrycode mktcode market geometry avg_price
## <chr> <dbl> <chr> <POINT [°]> <dbl>
## 1 AO 200 Luanda (13.23444 -8.838333) 1.80
## 2 BJ 401 Cotonou (2.433333 6.366667) 0.607
## 3 BJ 402 Malanville (3.393031 11.86195) 0.257
## 4 BJ 403 Natitingou (1.388616 10.31291) 0.459
## 5 BJ 404 Parakou (2.616667 9.35) 0.311
## 6 BW 110 Gaborone (25.90874 -24.65411) 0.2627 Calculate minimum distances
7.1 Coast line
# Load coastline data from the Natural Earth dataset as an sf object
coastline <- ne_coastline(scale = "medium", returnclass = "sf")# Calculate minimum distance from market places to coastline
merged_data$dist_to_coast <- apply(st_distance(merged_data, coastline), 1, min)
# Convert distance from meters to kilometers
merged_data$dist_to_coast <- merged_data$dist_to_coast / 10007.2 Check the minimum distance calculated
# Create graph to check if the minimum distance is correctly calculated
# Extract longitude and latitude from the geometry column of the sf object
merged_data <- merged_data %>%
mutate(longitude = st_coordinates(.)[,1],
latitude = st_coordinates(.)[,2])
# Filter the dataset to only include markets within Africa's bounding box
merged_data_africa <- merged_data %>%
filter(longitude > -25 & longitude < 60 & latitude > -40 & latitude < 40)
# Filter Coastline for Africa Only
coastline_africa <- coastline %>%
st_crop(xmin = -25, xmax = 60, ymin = -40, ymax = 40)
# Only include Africa
africa_map <- subset(world_map, long > -25 & long < 60 & lat > -40 & lat < 40)# Create map
ggplot() +
geom_polygon(data = africa_map, aes(x = long, y = lat, group = group),
fill = "gray90", color = "white") + # Africa background
geom_sf(data = coastline_africa, color = "blue", size = 0.5) + # Coastline in blue
geom_sf(data = merged_data_africa, aes(color = dist_to_coast), size = 2, alpha = 0.7) + # Market points
scale_color_viridis_c(option = "plasma") + # Better color scale
coord_sf(xlim = c(-25, 60), ylim = c(-40, 40), expand = FALSE) + # Restrict to Africa
theme_void() + # Removes grid lines and background
labs(title = "Market Distance to Coast", color = "Distance (km)") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) # Center title
7.3 Create scatter plot that relates average prices with coastline
# Create scatter plot that relates average prices with distance to coast line
ggplot(merged_data, aes(x = dist_to_coast, y = avg_price)) +
geom_point(alpha = 0.7, color = "blue") + # Scatter points
geom_smooth(method = "lm", color = "red", se = TRUE) + # Linear trend line
theme_minimal() +
labs(title = "Market Prices vs. Minimum Distance to Coast",
x = "Distance to Coast (km)",
y = "Average Market Price") +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))## `geom_smooth()` using formula = 'y ~ x'
The scatter plot does not reveal a clear pattern between market prices and distance to the coast. While the trend line suggests a slight positive relationship, the effect appears weak, with significant variation in prices at both short and long distances. The wide spread of data points and the broad confidence interval indicate that distance alone is not a strong predictor of market prices, suggesting that other factors, such as transportation infrastructure, local demand, or trade policies, may play a more significant role.
7.4 Roads
7.5 Calculate the minimum distance from markets to nearest roads
We now gather road data from Natural Earth. We filter on roads in Africa using the continent file and calculate the nearest road by first calculating the distance from each market to each road. Then we store the nearest roads with the shortest distances. We extract these roads and calculate the distance from market to nearest road in km. Finally we inspect our data to make sure the process ran smoothly.
# Define the URL to the shapefile on GitHub
shp_url <- "/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/Data/ne_10m_roads/ne_10m_roads.shp"
# Read the shapefile directly from the URL
road_data <- st_read(shp_url)## Reading layer `ne_10m_roads' from data source
## `/vsicurl/https://raw.githubusercontent.com/Viktoriag27/Geospatial_Assignment_2/main/Data/ne_10m_roads/ne_10m_roads.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 56600 features and 31 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -166.5325 ymin: -55.11212 xmax: 178.4191 ymax: 71.17768
## Geodetic CRS: WGS 84#Filter on raods in Africa
road_data<- road_data %>%
filter(continent == "Africa")
# Calculate the distance between each market point and all road geometries
distances <- st_distance(merged_data, road_data)
# Identify the nearest road for each market
nearest_road_index <- apply(distances, 1, which.min)
# Extract the nearest road geometries
nearest_roads <- road_data[nearest_road_index, ]
# Add the nearest road geometry to the merged_data dataframe
merged_data$nearest_road_geom <- nearest_roads$geometry
# Calculate the distance to the nearest road in km
merged_data$distance_to_nearest_road <- as.numeric(apply(distances, 1, min)/1000)
# View our changes in merged_data
head(merged_data)## Simple feature collection with 6 features and 8 fields
## Active geometry column: geometry
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1.388616 ymin: -24.65411 xmax: 25.90874 ymax: 11.86195
## Geodetic CRS: WGS 84
## # A tibble: 6 × 10
## ctrycode mktcode market geometry avg_price dist_to_coast
## <chr> <dbl> <chr> <POINT [°]> <dbl> <dbl>
## 1 AO 200 Luanda (13.23444 -8.838333) 1.80 1.57
## 2 BJ 401 Cotonou (2.433333 6.366667) 0.607 2.68
## 3 BJ 402 Malanville (3.393031 11.86195) 0.257 586.
## 4 BJ 403 Natitingou (1.388616 10.31291) 0.459 453.
## 5 BJ 404 Parakou (2.616667 9.35) 0.311 327.
## 6 BW 110 Gaborone (25.90874 -24.65411) 0.262 688.
## # ℹ 4 more variables: longitude <dbl>, latitude <dbl>,
## # nearest_road_geom <MULTILINESTRING [°]>, distance_to_nearest_road <dbl>For fun, we map the distance from markets to the nearest road.
# Create map
ggplot() +
geom_polygon(data = africa_map, aes(x = long, y = lat, group = group),
fill = "gray90", color = "white") + # Africa background
geom_sf(data = road_data, color = "green", size = 0.5) +
geom_sf(data = merged_data, aes(color = distance_to_nearest_road), size = 2, alpha = 0.7) + # Market points
scale_color_viridis_c(option = "cividis") + # Better color scale
coord_sf(xlim = c(-25, 60), ylim = c(-40, 40), expand = FALSE) + # Restrict to Africa
theme_void() + # Removes grid lines and background
labs(title = "Market Distance to the Nearest Road", color = "Distance (km)") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) # Center title
7.6 Create scatter plot that relates average prices with nearest road
# Create scatter plot that relates average prices with distance to nearest road
ggplot(merged_data, aes(x = distance_to_nearest_road, y = avg_price)) +
geom_point(alpha = 0.7, color = "blue") + # Scatter points
geom_smooth(method = "lm", color = "green", se = TRUE) + # Linear trend line
theme_minimal() +
labs(title = "Market Prices vs. Minimum Distance to Road",
x = "Distance to the Nearest Road (km)",
y = "Average Market Price") +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))## `geom_smooth()` using formula = 'y ~ x'
From this scatter plot, we detect no relationship between the average price of goods at the market and the distance to the nearest road. The trend line appears to have a positive slope but this is due to an outlier. Most markets are relatively close to roads and thus there are no identifiable trends in market pricing due to road access.
8 Airports
To gather data on airports in Africa, I used open source data from the afrimapr project and filtered on medium and large airports in Africa. I originally tried to use natural earth data; however, as there was no continent variable, I attempted to filter on airports in Africa by filtering on the bounding box. This created issues when I mapped my data as it also showed airports outside of Africa, such as on the Arabian Peninsula, that were within the bounding box. While it did not impact the distance calculation, it aesthetically did not look as organized or clear. I also tried to use st_within or st_intersect but encountered issues. Using the afrimapr data allowed for the simplest method and allowed me to map the markets and airports, which adds additional understanding to this project.
# Load the dataset for African airports
data("afriairports")
# We select relevant columns, filter on large and medium airports, ensure the geometry is treated as points, and that the correct CRS is being used
africa_airports <- afriairports %>%
select(name, country_name, geometry,type) %>%
filter(type %in% c("large_airport", "medium_airport")) %>%
st_as_sf() %>% # Ensure it remains an sf object
st_transform(crs = 4326) # Optional: Set CRS to WGS84 (if needed)
# View the first few rows to ensure its correct
print(head(africa_airports))## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -8.0363 ymin: -33.9648 xmax: 36.9278 ymax: 33.3675
## Geodetic CRS: WGS 84
## # A tibble: 6 × 4
## name country_name geometry type
## <chr> <chr> <POINT [°]> <chr>
## 1 OR Tambo International Airport South Africa (28.246 -26.1392) large_airpo…
## 2 Cairo International Airport Egypt (31.4056 30.1219) large_airpo…
## 3 Cape Town International Airport South Africa (18.6017 -33.9648) large_airpo…
## 4 Jomo Kenyatta International Airp… Kenya (36.9278 -1.31924) large_airpo…
## 5 Menara Airport Morocco (-8.0363 31.6069) medium_airp…
## 6 Mohammed V International Airport Morocco (-7.58997 33.3675) large_airpo…8.1 Calculate the minimum distance from markets to nearest airport
# Calculate minimum distance from market places to airport
merged_data$dist_to_airport <- apply(st_distance(merged_data, africa_airports), 1, min)
# Convert distance from meters to kilometers
merged_data$dist_to_airport <- merged_data$dist_to_airport / 1000
# View the results of our calculation
print(head(merged_data))## Simple feature collection with 6 features and 9 fields
## Active geometry column: geometry
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1.388616 ymin: -24.65411 xmax: 25.90874 ymax: 11.86195
## Geodetic CRS: WGS 84
## # A tibble: 6 × 11
## ctrycode mktcode market geometry avg_price dist_to_coast
## <chr> <dbl> <chr> <POINT [°]> <dbl> <dbl>
## 1 AO 200 Luanda (13.23444 -8.838333) 1.80 1.57
## 2 BJ 401 Cotonou (2.433333 6.366667) 0.607 2.68
## 3 BJ 402 Malanville (3.393031 11.86195) 0.257 586.
## 4 BJ 403 Natitingou (1.388616 10.31291) 0.459 453.
## 5 BJ 404 Parakou (2.616667 9.35) 0.311 327.
## 6 BW 110 Gaborone (25.90874 -24.65411) 0.262 688.
## # ℹ 5 more variables: longitude <dbl>, latitude <dbl>,
## # nearest_road_geom <MULTILINESTRING [°]>, distance_to_nearest_road <dbl>,
## # dist_to_airport <dbl>We map the distance from markets to the nearest airport.
library(ggplot2)
# Create map
ggplot() +
geom_polygon(data = africa_map, aes(x = long, y = lat, group = group),
fill = "gray90", color = "white") + # Africa background
geom_sf(data = africa_airports, color = "red", size = 0.5) + # Airports in blue
geom_sf(data = merged_data, aes(color = dist_to_airport), size = 2, alpha = 0.7) + # Market points
scale_color_viridis_c(option = "magma") + # Better color scale
coord_sf(xlim = c(-25, 60), ylim = c(-40, 40), expand = FALSE) + # Restrict to Africa
theme_void() + # Removes grid lines and background
labs(title = "Market Distance to Airport", color = "Distance (km)") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) # Center title
8.2 Create scatter plot that relates average prices with nearest airport
# Create scatter plot that relates average prices with distance to airports
ggplot(merged_data, aes(x = dist_to_airport, y = avg_price)) +
geom_point(alpha = 0.7, color = "blue") + # Scatter points
geom_smooth(method = "lm", color = "pink", se = TRUE) + # Linear trend line
theme_minimal() +
labs(title = "Market Prices vs. Minimum Distance to Airport",
x = "Distance to Airport (km)",
y = "Average Market Price") +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))## `geom_smooth()` using formula = 'y ~ x'
Observing our results from the scatter plot, it is unclear if there is a relationship between the average price of goods at the market and the distance to the nearest airport. The trend line moves slightly downward, indicating a slight negative correlation as with price as distance increases. However, our data contains many higher and lower contradictory outliers.